%% 
% 'a' needs to be Nx2 array, where first column is temperature, second is
% growthrate
load('gr_data_all_20230413_v2.mat');
close all
a = gr_data_all;
T = unique(a(:,1));

G = [];
Gerr = [];
for k=1:6
    g = reshape(a(:,k+1),[length(a)/length(T) length(T)]);
    gmean = nanmean(g,1)';
    gstd = nanstd(g,1)';

    G = [G gmean];
    Gerr = [Gerr gstd];

end


%% Plot data
T1plot = 16;
T2plot = 50;

% clrs = min(distinguishable_colors(6)*3, 1);
clrs = (cbrewer('qual', 'Set1', 6));
clrs = (cbrewer('qual', 'Paired',6));
strainnames = {'MG1655', 'BW25113', 'CS109', 'BL21', 'dnaK', 'uup'};

ms = 6;
figure;

toplot = [1:4];

for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
    errorbar(T,g_plot,gerr_plot,'ko', 'markersize', ms, 'markerfacecolor', clrs(k,:));
    hold on;
    
end

set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('Temperature (°C)');
ylabel('Growth rate (1/h)');
box off;
legend(strainnames{toplot})
xlim([T1plot T2plot])

%% Plot Arrhenius data

Tmin = 25;
Tmax = 37;
clc;
Ea_tot = [];
Ea_err_tot = [];


figure;
for k=toplot
    g_plot = G(:,k);
    gerr_plot = Gerr(:,k);
    errorbar(1./(T+273),log(g_plot), gerr_plot./g_plot, 'ko', 'markersize', ms, 'markerfacecolor', clrs(k,:));
    hold on;

    Tfit = T(T<=Tmax & T>=Tmin);
    gfit = g_plot(T<=Tmax & T>=Tmin);
    gfit_err = gerr_plot(T<=Tmax & T>=Tmin);
    
    Tfit = Tfit(~isnan(gfit));
    gfit_err = gfit_err(~isnan(gfit));
    gfit = gfit(~isnan(gfit));
    

    % Linear fit
    [f1,f2] = fit(1./(Tfit + 273.15), log(gfit), 'poly1');

    % Weighted linear fit
    ft = fittype('poly1');
    cf = fit(1./(Tfit + 273.15), log(gfit),ft,'Weight',gfit_err./gfit);
    
    % Weighted linear fit using fitnlm
    Arrhenius_fit = @(b,x) b(1) + b(2).*x;
    x0 = [30 -8000];
    f = fitnlm(1./(Tfit + 273.15), log(gfit),Arrhenius_fit,x0, 'Weights', gfit_err./gfit);
    hold on;
    plot(1./(Tfit + 273.15), f.feval(1./(Tfit + 273.15)), 'color', clrs(k,:), 'linewidth', 2)

    Ea_tot = [Ea_tot; -f.Coefficients.Estimate(2)/298/1.688];
    Ea_err_tot = [Ea_err_tot; f.Coefficients.SE(2)/298/1.688];
end
set(gcf, "Position", [0 0 400 300]);
set(gca, 'FontSize', 20);
xlabel('1000/T (1/K)');
ylabel('Growth rate (1/h)');
xlim([1/(T2plot + 273.15) 1/(T1plot + 273.15)])
box off;


%% Plot bars with errors of activation energies (in units of kcal/mol)
x = 1:length(strainnames);
figure;

hold on
for k = 1:length(Ea_tot)
    er = errorbar(x(k),Ea_tot(k),0,Ea_err_tot(k), 'markeredgecolor', clrs(toplot(k),:)); 
    bar(k,Ea_tot(k), 'facecolor', clrs(toplot(k),:));
    hold on;
    er.Color = [0 0 0];                            

end

box off;
set(gcf, 'Position', [0 0 400 150]);
set(gca,'fontsize', 20);
ylabel('Activation energy');
ylim([0 20])